par(mfrow=c(2, 2))
data1_20 <- read.csv("data/时间序列分析——基于R（第2版）案例数据/csv/A1_20.csv")
ts_x <- ts(data1_20$x, frequency = 4, start=c(1962, 1))

# 时序图
plot(ts_x, type="o", pch=5, col='#39CBB4')

# 1阶4步差分
diff_x <- diff(diff(ts_x, 4))
plot(diff_x)

# ADF检验
# install.packages("aTSA")
library(aTSA)
adf.test(diff_x, nlag=3)

# 白噪声，纯随机检验（<0.05为非白噪声）
for( k in 1:3) print(Box.test(diff_x, lag=6*k, type="Ljung-Box"))

# 自相关图
acf(diff_x, lag.max=30)
pacf(diff_x, lag.max=30)

# 拟合加法ARIMA模型
fit <- arima(ts_x, order=c(4, 1, 0), seasonal = list(order=c(0, 1, 0), period=4), transform.pars = F, fixed=c(NA, 0, 0, NA))
fit

# 模型显著性检验
ts.diag(fit)

# 三年预测
# install.packages("forecast")
library(forecast)
x_fore <- forecast(fit, h=12)
x_fore

# 绘制效果图
plot(x_fore, lty=2)
lines(fitted(x_fore), col=4)

